unify access to ellipsoid parameters (#1312)

* define constants for often used spherioid parms.

* add function to calculate semi-minor axis.

* use semi major axis function more.

* have semi_minor_axis consume an ellipse

* compute semi minor axis in GPS_Ellipse class.

add constants for other ellipses that are used directly.

diff --git a/jeeps/gpsdatum.h b/jeeps/gpsdatum.h
index 2cfbb366aae0d85cc0b87034082d230c3c311cbf..eea7a361200c96c3097a1070a4da71a1fa9baf8b 100644 (file)
--- a/jeeps/gpsdatum.h
+++ b/jeeps/gpsdatum.h
@@ -6,14 +6,33 @@ struct GPS_Ellipse {
   const char*   name;
   double a;
   double invf;
+
+  constexpr double b() const {
+    return a - a/invf;
+  }
 };
 
+// EPSG:7001
+constexpr GPS_Ellipse Airy1830_Ellipse = { "Airy 1830", 6377563.396, 299.3249646 };
+
+// EPSG:7002
+constexpr GPS_Ellipse Airy1830Modified_Ellipse = { "Airy 1830 Modified", 6377340.189, 299.3249646 };
+
+// EPSG:7004
+constexpr GPS_Ellipse Bessel1841_Ellipse = { "Bessel 1841", 6377397.155, 299.1528128 };
+
+// EPSG:7019
+constexpr GPS_Ellipse GRS80_Ellipse = { "GRS80", 6378137.000, 298.257222101 };
+
+// EPSG:4326
+constexpr GPS_Ellipse WGS84_Ellipse = { "WGS84", 6378137.000, 298.257223563 };
+
 const GPS_Ellipse GPS_Ellipses[]= {
-  { "Airy 1830",               6377563.396, 299.3249646 },
-  { "Airy 1830 Modified",      6377340.189, 299.3249646 },
+  Airy1830_Ellipse,
+  Airy1830Modified_Ellipse,
   { "Australian National",     6378160.000, 298.25 },
   { "Bessel 1841 (Namibia)",   6377483.865, 299.1528128 },
-  { "Bessel 1841",             6377397.155, 299.1528128 },
+  Bessel1841_Ellipse,
   { "Clarke 1866",             6378206.400, 294.9786982 },
   { "Clarke 1880",             6378249.145, 293.465 },
   { "Everest (India 1830)",    6377276.345, 300.8017 },
@@ -30,12 +49,12 @@ const GPS_Ellipse GPS_Ellipses[]= {
   { "Krassovsky 1940",         6378245.000, 298.3 },
   { "GRS67",                   6378160.000, 6356774.516 },
   { "GRS75",                   6378140.000, 6356755.288 },
-  { "GRS80",                   6378137.000, 298.257222101 },
+  GRS80_Ellipse,
   { "S. American 1969",        6378160.000, 298.25 },
   { "WGS60",                   6378165.000, 298.3 },
   { "WGS66",                   6378145.000, 298.25 },
   { "WGS72",                   6378135.000, 298.26 },
-  { "WGS84",                   6378137.000, 298.257223563 },
+  WGS84_Ellipse,
   { "Clarke 1880 (Benoit)",    6378300.789, 293.466 },
 };
 

diff --git a/jeeps/gpsmath.cc b/jeeps/gpsmath.cc
index c97b6ad48c21aadcfb315c458085b9becaa02bae..6edfb5a145d0cdf2636f95d2a555d268169d478c 100644 (file)
--- a/jeeps/gpsmath.cc
+++ b/jeeps/gpsmath.cc
@@ -429,8 +429,8 @@ void GPS_Math_XYZ_To_LatLonH(double* phi, double* lambda, double* H,
 void GPS_Math_Airy1830LatLonH_To_XYZ(double phi, double lambda, double H,
                                      double* x, double* y, double* z)
 {
-  double a = 6377563.396;
-  double b = 6356256.910;
+  constexpr double a = Airy1830_Ellipse.a;
+  constexpr double b = Airy1830_Ellipse.b();
 
   GPS_Math_LatLonH_To_XYZ(phi,lambda,H,x,y,z,a,b);
 
@@ -456,8 +456,8 @@ void GPS_Math_Airy1830LatLonH_To_XYZ(double phi, double lambda, double H,
 void GPS_Math_WGS84LatLonH_To_XYZ(double phi, double lambda, double H,
                                   double* x, double* y, double* z)
 {
-  double a = 6378137.000;
-  double b = 6356752.314245;
+  constexpr double a = WGS84_Ellipse.a;
+  constexpr double b = WGS84_Ellipse.b();
 
   GPS_Math_LatLonH_To_XYZ(phi,lambda,H,x,y,z,a,b);
 
@@ -483,8 +483,8 @@ void GPS_Math_WGS84LatLonH_To_XYZ(double phi, double lambda, double H,
 void GPS_Math_XYZ_To_Airy1830LatLonH(double* phi, double* lambda, double* H,
                                      double x, double y, double z)
 {
-  double a = 6377563.396;
-  double b = 6356256.910;
+  constexpr double a = Airy1830_Ellipse.a;
+  constexpr double b = Airy1830_Ellipse.b();
 
   GPS_Math_XYZ_To_LatLonH(phi,lambda,H,x,y,z,a,b);
 
@@ -509,8 +509,8 @@ void GPS_Math_XYZ_To_Airy1830LatLonH(double* phi, double* lambda, double* H,
 void GPS_Math_XYZ_To_WGS84LatLonH(double* phi, double* lambda, double* H,
                                   double x, double y, double z)
 {
-  double a = 6378137.000;
-  double b = 6356752.314245;
+  constexpr double a = WGS84_Ellipse.a;
+  constexpr double b = WGS84_Ellipse.b();
 
   GPS_Math_XYZ_To_LatLonH(phi,lambda,H,x,y,z,a,b);
 
@@ -643,8 +643,8 @@ void GPS_Math_Airy1830M_LatLonToINGEN(double phi, double lambda, double* E,
   double F0      = 1.000035;
   double phi0    = 53.5;
   double lambda0 = -8.;
-  double a       = 6377340.189;
-  double b       = 6356034.447;
+  constexpr double a       = Airy1830Modified_Ellipse.a;
+  constexpr double b       = Airy1830Modified_Ellipse.b();
 
   GPS_Math_LatLon_To_EN(E,N,phi,lambda,N0,E0,phi0,lambda0,F0,a,b);
 
@@ -674,8 +674,8 @@ void GPS_Math_Airy1830LatLonToNGEN(double phi, double lambda, double* E,
   double F0      = 0.9996012717;
   double phi0    = 49.;
   double lambda0 = -2.;
-  double a       = 6377563.396;
-  double b       = 6356256.910;
+  constexpr double a       = Airy1830_Ellipse.a;
+  constexpr double b       = Airy1830_Ellipse.b();
 
   GPS_Math_LatLon_To_EN(E,N,phi,lambda,N0,E0,phi0,lambda0,F0,a,b);
 
@@ -704,7 +704,7 @@ int32_t GPS_Math_WGS84_To_Swiss_EN(double lat, double lon, double* E,
   const double lambda0 = 7.43958333;
   const double E0 = 600000.0;
   const double N0 = 200000.0;
-  double phi, lambda, alt, a, b;
+  double phi, lambda, alt;
 
   if (lat < 44.89022757) {
     return 0;
@@ -713,9 +713,8 @@ int32_t GPS_Math_WGS84_To_Swiss_EN(double lat, double lon, double* E,
     return 0;
   }
 
-  assert(strcmp(GPS_Ellipses[4].name, "Bessel 1841") == 0);
-  a = GPS_Ellipses[4].a;
-  b = a - (a / GPS_Ellipses[4].invf);
+  constexpr double a = Bessel1841_Ellipse.a;
+  constexpr double b = Bessel1841_Ellipse.b();
 
   GPS_Math_WGS84_To_Known_Datum_M(lat, lon, 0, &phi, &lambda, &alt, 123);
   GPS_Math_Swiss_LatLon_To_EN(phi, lambda, E, N, phi0, lambda0, E0, N0, a, b);
@@ -742,11 +741,10 @@ void GPS_Math_Swiss_EN_To_WGS84(double E, double N, double* lat, double* lon)
   const double lambda0 = 7.43958333;
   const double E0 = 600000.0;
   const double N0 = 200000.0;
-  double phi, lambda, alt, a, b;
+  double phi, lambda, alt;
 
-  assert(strcmp(GPS_Ellipses[4].name, "Bessel 1841") == 0);
-  a = GPS_Ellipses[4].a;
-  b = a - (a / GPS_Ellipses[4].invf);
+  constexpr double a = Bessel1841_Ellipse.a;
+  constexpr double b = Bessel1841_Ellipse.b();
 
   GPS_Math_Swiss_EN_To_LatLon(E, N, &phi, &lambda, phi0, lambda0, E0, N0, a, b);
   GPS_Math_Known_Datum_To_WGS84_M(phi, lambda, 0, lat, lon, &alt, 123);
@@ -1112,7 +1110,7 @@ int32_t GPS_Math_WGS84_To_ICS_EN(double lat, double lon, double* E,
   int32_t ellipse = GPS_Datums[datum].ellipse;
 
   a = GPS_Ellipses[ellipse].a;
-  b = a - (a / GPS_Ellipses[ellipse].invf);
+  b = GPS_Ellipses[ellipse].b();
 
   GPS_Math_WGS84_To_Known_Datum_M(lat, lon, 0, &phi, &lambda, &alt, datum);
   GPS_Math_Cassini_LatLon_To_EN(phi, lambda, E, N,
@@ -1148,7 +1146,7 @@ void GPS_Math_ICS_EN_To_WGS84(double E, double N, double* lat, double* lon)
   int32_t ellipse = GPS_Datums[datum].ellipse;
 
   a = GPS_Ellipses[ellipse].a;
-  b = a - (a / GPS_Ellipses[ellipse].invf);
+  b = GPS_Ellipses[ellipse].b();
 
   GPS_Math_Cassini_EN_To_LatLon(E, N, &phi, &lambda, phi0, lambda0,
                                 E0, N0, a, b);
@@ -1307,8 +1305,8 @@ void GPS_Math_NGENToAiry1830LatLon(double E, double N, double* phi,
   double F0      = 0.9996012717;
   double phi0    = 49.;
   double lambda0 = -2.;
-  double a       = 6377563.396;
-  double b       = 6356256.910;
+  constexpr double a = Airy1830_Ellipse.a;
+  constexpr double b = Airy1830_Ellipse.b();
 
   GPS_Math_EN_To_LatLon(E,N,phi,lambda,N0,E0,phi0,lambda0,F0,a,b);
 
@@ -1337,8 +1335,8 @@ void GPS_Math_INGENToAiry1830MLatLon(double E, double N, double* phi,
   double F0      = 1.000035;
   double phi0    = 53.5;
   double lambda0 = -8.;
-  double a       = 6377340.189;
-  double b       = 6356034.447;
+  constexpr double a = Airy1830Modified_Ellipse.a;
+  constexpr double b = Airy1830Modified_Ellipse.b();
 
   GPS_Math_EN_To_LatLon(E,N,phi,lambda,N0,E0,phi0,lambda0,F0,a,b);
 
@@ -1539,15 +1537,13 @@ void GPS_Math_Known_Datum_To_WGS84_M(double Sphi, double Slam, double SH,
 {
   double Sa;
   double Sif;
-  double Da;
-  double Dif;
   double x;
   double y;
   double z;
   int32_t idx;
 
-  Da  = 6378137.0;
-  Dif = 298.257223563;
+  constexpr double Da  = WGS84_Ellipse.a;
+  constexpr double Dif = WGS84_Ellipse.invf;
 
   idx  = GPS_Datums[n].ellipse;
   Sa   = GPS_Ellipses[idx].a;
@@ -1581,8 +1577,6 @@ void GPS_Math_WGS84_To_Known_Datum_M(double Sphi, double Slam, double SH,
                                      double* Dphi, double* Dlam, double* DH,
                                      int32_t n)
 {
-  double Sa;
-  double Sif;
   double Da;
   double Dif;
   double x;
@@ -1590,8 +1584,8 @@ void GPS_Math_WGS84_To_Known_Datum_M(double Sphi, double Slam, double SH,
   double z;
   int32_t idx;
 
-  Sa  = 6378137.0;
-  Sif = 298.257223563;
+  constexpr double Sa  = WGS84_Ellipse.a;
+  constexpr double Sif = WGS84_Ellipse.invf;
 
   idx  = GPS_Datums[n].ellipse;
   Da   = GPS_Ellipses[idx].a;
@@ -1628,9 +1622,6 @@ void GPS_Math_Known_Datum_To_WGS84_C(double Sphi, double Slam, double SH,
   double Sa;
   double Sif;
   double Sb;
-  double Da;
-  double Dif;
-  double Db;
   double x;
   double y;
   double z;
@@ -1639,9 +1630,8 @@ void GPS_Math_Known_Datum_To_WGS84_C(double Sphi, double Slam, double SH,
   double sy;
   double sz;
 
-  Da  = 6378137.0;
-  Dif = 298.257223563;
-  Db  = Da - (Da / Dif);
+  constexpr double Da  = WGS84_Ellipse.a;
+  constexpr double Db  = WGS84_Ellipse.b();
 
   idx  = GPS_Datums[n].ellipse;
   Sa   = GPS_Ellipses[idx].a;
@@ -1682,23 +1672,19 @@ void GPS_Math_WGS84_To_Known_Datum_C(double Sphi, double Slam, double SH,
                                      double* Dphi, double* Dlam, double* DH,
                                      int32_t n)
 {
-  double Sa;
-  double Sif;
   double Da;
   double Dif;
   double x;
   double y;
   double z;
   int32_t idx;
-  double Sb;
   double Db;
   double dx;
   double dy;
   double dz;
 
-  Sa  = 6378137.0;
-  Sif = 298.257223563;
-  Sb   = Sa - (Sa / Sif);
+  constexpr double Sa  = WGS84_Ellipse.a;
+  constexpr double Sb  = WGS84_Ellipse.b();
 
   idx  = GPS_Datums[n].ellipse;
   Da   = GPS_Ellipses[idx].a;
@@ -1803,9 +1789,7 @@ void GPS_Math_Known_Datum_To_Known_Datum_C(double Sphi, double Slam, double SH,
     double* DH, int32_t n1, int32_t n2)
 {
   double Sa;
-  double Sif;
   double Da;
-  double Dif;
   double x1;
   double y1;
   double z1;
@@ -1824,8 +1808,7 @@ void GPS_Math_Known_Datum_To_Known_Datum_C(double Sphi, double Slam, double SH,
 
   idx1  = GPS_Datums[n1].ellipse;
   Sa    = GPS_Ellipses[idx1].a;
-  Sif   = GPS_Ellipses[idx1].invf;
-  Sb    = Sa - (Sa / Sif);
+  Sb    = GPS_Ellipses[idx1].b();
 
   x1    = GPS_Datums[n1].dx;
   y1    = GPS_Datums[n1].dy;
@@ -1833,8 +1816,7 @@ void GPS_Math_Known_Datum_To_Known_Datum_C(double Sphi, double Slam, double SH,
 
   idx2  = GPS_Datums[n2].ellipse;
   Da    = GPS_Ellipses[idx2].a;
-  Dif   = GPS_Ellipses[idx2].invf;
-  Db    = Da - (Da / Dif);
+  Db    = GPS_Ellipses[idx2].b();
 
   x2    = GPS_Datums[n2].dx;
   y2    = GPS_Datums[n2].dy;
@@ -2116,8 +2098,6 @@ int32_t GPS_Math_NAD83_To_UTM_EN(double lat, double lon, double* E,
   double N0;
   double E0;
   double F0;
-  double a;
-  double b;
 
   if (!GPS_Math_LatLon_To_UTM_Param(lat,lon,zone,zc,&lambda0,&E0,
                                     &N0,&F0)) {
@@ -2126,9 +2106,8 @@ int32_t GPS_Math_NAD83_To_UTM_EN(double lat, double lon, double* E,
 
   phi0 = 0.0;
 
-  assert(strcmp(GPS_Ellipses[21].name, "GRS80") == 0);
-  a = GPS_Ellipses[21].a;
-  b = a - (a / GPS_Ellipses[21].invf);
+  constexpr double a = GRS80_Ellipse.a;
+  constexpr double b = GRS80_Ellipse.b();
 
   GPS_Math_LatLon_To_EN(E,N,lat,lon,N0,E0,phi0,lambda0,F0,a,b);
 
@@ -2294,7 +2273,7 @@ int32_t GPS_Math_Known_Datum_To_UTM_EN(double lat, double lon, double* E,
 
   idx  = GPS_Datums[n].ellipse;
   a = GPS_Ellipses[idx].a;
-  b = a - (a / GPS_Ellipses[idx].invf);
+  b = GPS_Ellipses[idx].b();
 
   GPS_Math_LatLon_To_EN(E,N,lat,lon,N0,E0,phi0,lambda0,F0,a,b);
 
@@ -2517,7 +2496,7 @@ void GPS_Math_UTM_EN_to_LatLon(int ReferenceEllipsoid,
 //found at http://www.gpsy.com/gpsinfo/geotoutm/index.html
 
   double k0 = 0.9996;
-  double a, b;
+  double a, f;
   double eccSquared;
   double eccPrimeSquared;
   double e1;
@@ -2526,8 +2505,8 @@ void GPS_Math_UTM_EN_to_LatLon(int ReferenceEllipsoid,
   double x, y;
 
   a = GPS_Ellipses[ReferenceEllipsoid].a;
-  b = 1 / GPS_Ellipses[ReferenceEllipsoid].invf;
-  eccSquared = b * (2.0 - b);
+  f = 1 / GPS_Ellipses[ReferenceEllipsoid].invf;
+  eccSquared = f * (2.0 - f);
   e1 = (1-sqrt(1-eccSquared))/(1+sqrt(1-eccSquared));
 
   x = UTMEasting - E0; //remove false easting

diff --git a/reference/grid-bng~csv.gpx b/reference/grid-bng~csv.gpx
index 0c361caec8ecb1379295a34b853dfb24d0ee61a7..ddcb600ea1aeea8e8d6fb56ebccbcdfaa28700e7 100644 (file)
--- a/reference/grid-bng~csv.gpx
+++ b/reference/grid-bng~csv.gpx
@@ -1,8 +1,8 @@
 <?xml version="1.0" encoding="UTF-8"?>
 <gpx version="1.0" creator="GPSBabel - https://www.gpsbabel.org" xmlns="http://www.topografix.com/GPX/1/0">
   <time>1970-01-01T00:00:00Z</time>
-  <bounds minlat="50.207517637" minlon="-7.489235773" maxlat="60.150151215" maxlon="1.296391371"/>
-  <wpt lat="58.007809302" lon="-6.751777638">
+  <bounds minlat="50.207517638" minlon="-7.489235772" maxlat="60.150151215" maxlon="1.296391371"/>
+  <wpt lat="58.007809302" lon="-6.751777637">
     <time>2008-08-23T20:56:12Z</time>
     <name>AlineLodge</name>
     <cmt>Aline Lodge</cmt>
@@ -16,7 +16,7 @@
     <desc>Caernarfon</desc>
     <sym>City (Small)</sym>
   </wpt>
-  <wpt lat="50.207517637" lon="-5.295409028">
+  <wpt lat="50.207517638" lon="-5.295409027">
     <time>2008-08-23T20:55:39Z</time>
     <name>Camborne</name>
     <cmt>Camborne</cmt>
@@ -37,28 +37,28 @@
     <desc>Forfar</desc>
     <sym>City (Small)</sym>
   </wpt>
-  <wpt lat="55.426025760" lon="-2.790531785">
+  <wpt lat="55.426025761" lon="-2.790531785">
     <time>2008-08-23T20:55:25Z</time>
     <name>Hawick</name>
     <cmt>Hawick</cmt>
     <desc>Hawick</desc>
     <sym>City (Small)</sym>
   </wpt>
-  <wpt lat="57.621914847" lon="-7.489235773">
+  <wpt lat="57.621914847" lon="-7.489235772">
     <time>2008-08-23T20:56:17Z</time>
     <name>Hosta</name>
     <cmt>Hosta</cmt>
     <desc>Hosta</desc>
     <sym>City (Small)</sym>
   </wpt>
-  <wpt lat="57.007373614" lon="-6.271139684">
+  <wpt lat="57.007373614" lon="-6.271139683">
     <time>2008-08-23T20:56:24Z</time>
     <name>IsleofRhum</name>
     <cmt>Isle of Rhum</cmt>
     <desc>Isle of Rhum</desc>
     <sym>City (Small)</sym>
   </wpt>
-  <wpt lat="60.150151215" lon="-1.142586362">
+  <wpt lat="60.150151215" lon="-1.142586363">
     <time>2008-08-23T20:59:39Z</time>
     <name>Lerwick</name>
     <cmt>Lerwick</cmt>